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Hillel Tal-Ezer 

School of Mathematical Sciences 
Tel-Aviv University, Tel-Aviv, Israel 


Abstract 

A pseudospectral numerical scheme for solving linear, periodic, hyperbolic 
problems is described. It has infinite accuracy both in time and in space. 
The high accuracy in time is achieved without increasing the computational 
work and memory space which is needed for a regular, one step explicit 
scheme. The algorithm is shown to be optimal in the sense that among all the 
explicit algorithms of a certain class it requires the least amount of work to 
achieve a certain given resolution. The class of algorithms referred to 
consists of all explicit schemes which may be represented as a polynomial in 
the spatial operator. 
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1. Introduction 


In recent years, it has been shown that spectral methods can provide a 
very useful tool for the solution of time dependent partial differential 
equations. A standard scheme uses spectral methods to approximate the space 
derivatives and a finite difference approach to march the solution in time. 
This tactic results in an unbalanced scheme; it has infinite accuracy in space 
and finite accuracy in time. It is obvious that the overall accuracy is 
influenced strongly by the relatively poor approximation of the time 
derivative. 

In this paper we present an alternative approach that also yields spectral 
accuracy in time and is optimal in terms of efficiency. 

The finite difference approach for the time discretization of the P.D.E. 
is discussed in Section 2. In Section 3 we present the new approach for 
marching the solution in time in order to get an overall infinite accuracy. 

The method presented in Section 3 is based on expanding the evolution 
operator by orthogonal polynomials. In Sections 4 and 5 we discuss the 
resolution and stability properties of the method. The scheme is compared to 
the leap-frog type schemes in order to clarify its properties. In Section 6 
we give a proof of the infinite accuracy of our approach. Section 7 presents 
the algorithm in detail and in Section 8 we demonstrate its validity for the 
variable coefficients case. Section 9 concludes the discussion by presenting 
numerical results which confirm the theoretical results developed in the 
previous sections. 
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2. Finite Difference Approach 

Consider the differential equation 


( 2 . 1 ) 


U - GU = 0 0 < x < 2tt 

U(x,0) = U°(x) 


where G is a linear spatial differential operator. We assume that the 
coefficients of the derivatives appearing in G are time independent and 
2 it - periodic. Suppose further that (2.1) is discretized in space by using 
the pseudospectral Fourier method, [2], [4], This involves seeking a 

trigonometric polynomial Ujj(x) of degree N that satisfies 


( 2 . 2 ) 


3U 

Tt 


N 


Wn 


0 


U N (x,0) = P N U°(x) 


where for any function f(x), P N f(x) is its trigonometric interpolant at the 
collocation points 


more precisely, 
(2.3) 


V* 


j - 0,1 ,• • • ,2N-1; 


N 


P f(x) = l a e 
N k=-N 


Ikx 


, 2N-1 

= 2? I £ W 


-ikx. 

e 3 


j=0 


where 


3 


The solution of (2.2) is given by 

< 2 * 4 > V X,t) = ex P( tP N GP ^ U ° (x) * 

Except for very simple operators G, it is impractical to construct the 

exponential matrix exp( tP^GP^) explicity. Usually an approximation to the 

exponent is used. Most frequently an explicit or implicit finite difference 

scheme is used to march the solution over a time step At. All these 

algorithms are based on a Taylor expansion of the evolution operator 

exp( tP^GP^) . Essentially, the scalar function e z is expanded either as a 

m 

Taylor series of the for £ z fZ\ or by a Pade approximation 

Z = 0 


(2.5) 


z 


e 


P 


l b 
£=0 


z 


z 

z 



z=o 


z 


z 

z 


where b^, c^, are so chosen that the expansion of the right side of (2.5) 
agrees with the Taylor expansion and then z is replaced by the matrix 
4t(P„G p N ). 

For example, in the case of the modified Euler scheme one advances from 
the time leve n to n + 1 by using 


1 + AtP N GP N 


+ 



2 


to approximate exp[At(P N GP N ) ] . Let V n be the approximation to U N (n*At). 
Then 
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( 2 . 6 ) 

or 

(2.7) 
where 

( 2 . 8 ) 


V n+1 = [l+At(P N GP N )+^(P N GP N ) 2 ]v n 
V n = (l+AtG N + ^|-G 2 ) n U N (0) 

g n = p n gp n* 


Equation (2.7) can be rewritten as 


(2.9) 


v "- 1 


or 


y n 


2n 


v “ = 1 J 0 P k^V) |V 0) 


where 


- p k 

e k = -k 
n 


These types of approximations result in a limitation on the allowable time 
step At since Taylor expansion possesses high accuracy for small At and 
this accuracy deteriorates rapidly as At increases. 

For problems in which the solution changes over a time scale which is 
comparable to the spatial scale, as in hyperbolic problems, the limitation on 
At can make the scheme impractical. We therefore explore other possible 
expansions that do not suffer from this drawback. A natural candidate is an 
expansion based on orthogonal polynomials. We thus restrict our discussion to 
polynomial schemes, i.e. , to schemes that employ an algorithm which 
approximates the numerical solution at time level t in the following way 
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V(x,t) = l a (G t) k V°(x) 
k=0 N 

The modified Euler algorithm (2.7) is an example of such a scheme, as well as 
are most explicit methods. 


3. Orthogonal Polynomial Approach 

We start by explaining how the new method is constructed in the case of 
the simple hyperbolic equation 

U - aU =0 
t x 

(3.1) 0 <C x _< 3 TT 

U(x,0) = Uq(x) 

U(2tt ,t) = U(0,t) t > 0 


where a is constant. 

The semi -discrete pseudospectral Fourier method can be written in the form 


(3.2) 


3 V 
3 1 


= g n v 


V(t=0) = v° 


where V(t) is the column vector 


( v ( x 0 )> v ( x i) »••• » v (x 2n _ 1 )) T 


and 
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v° - (O 0 (x 0 )- — *’ , 0 ( X 2N-1»' 


The 2N x 2N matrix G N is given by 


(3.3) 



(fC-D^ctgiU 

I • 


j * k 
j = k . 


(In practice G^V is calculated using two fast Fourier transforms.) Gjg is a 
skew-symmetric matrix and therefore has a complete set of 2N eigenvectors 
which will be deonted by k = 1,2, •••2N. Let 


2N 

- l b 

k=l 


k“k 


then 

(3.4) 


G KT t 2N G„t 2N X, t 

N-rjO V , N— Vv k — 

e v = l b, e t*. = l b, e w, 

k=l K K k=l 


V 


where X^ are the eigenvalues of G ^ corresponding to In our case X 

are purely imaginery. Let be a polynomial approximation of 

exp(G N t) of the form 



then 




2N _ 

JA h JvK 

k=l 


and therefore 


ll[e V -»JV)^N 2 - ri\| 2 |>\(V )| 2 

k=l 
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if are arbitrary (see however the remark at the end of this section), 

(3.5) | le^ -H m (G N t)|| 2 = max|e k | 2 £ max|e z -H m (z) | 2 

k z 

where 

z € [iat(N-l ) , iat(N-l) ] 

We therefore seek a polynomial approximation with real coefficients to the 
function e z that will minimize the expression on the R.H.S. of (3.5). 

Define 

(3.6) R = |at(N-l) | 

(3.7) 0 = -iz/R ( |0 | < 1) 

Then 

l eZ "V z) | 2 = |e i0R -H m (i9R)| 2 = 

|cos(0R)-H R (0R) I 2 + lsin(0R)-H I (0R) I 2 
m 1 m 1 

where 

(3.8) H (i0R) = H R (0R) + iH I (0R) 

m m m 

R X 

(Hm, are polynomials with real coefficients.) 

The polynomials that minimize (3.8) are the "best approximation" to cos(9R) 
and sin(0R). It is known that Chebyshev polynomials provide an approximation 
which is "almost" as good as the "best approximation". In fact, we can quote 
the following result [7]. 
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Theorem : Suppose that f 6 c[-l,l] and S^(f) = llf-q^l^ where 


q n is 


the least-square approximation with respect to the weight function ( l-z^) ^ 


then 


S n (f) < (4+ -^log n)E n (f) 


ir 


(E^Cf) = q n is the best polynomial approximation.) 

It follows that the improved accuracy of the best polynomial approximation 
does not make up for the added computational complexity. Taking H^CBR), 
H^(0R) as the Chebyshev polynomial approximations to the trigonometric 
functions we have [8] 


(3.9) 


CO 

H*(0R) = J (R) +2 l (-l) k J (R)T (0) 
m o , - Zk Zk 

k=l 


- 2 i C-l> k J 2t+1 «W0) 

k=l 


where J^CR) is Vessel function of order K. Hence 


m , 

(3.10) H m (i9R) = ^ (i) c k J fc (R)T k (0) 

k=0 


c 0 = !> c k = 2 k 1 1 

Since (3.7) we have 

(3.11) T k (0) = T k^“i w ^ w = Z / R w € t“i>i] 
Define 

(3.12) Q k (w) = (i) k T k (-iw) ; 


using the recurrence relation satisfied by Chebyshev polynomials 



9 


(3 * 13) T k+l (x) = 2xT k-l (x) " T k-l (x) T 0 (x) = V*) - x 

it is easily verified that Q^(w) satisfies the following recurrence relation 


(3.14) 


Qk+i ( w ) = 2wQ k (w) + Q^jCw) 
Q q (w) =1, QjCw) = w 


Thus, Q k 's are polynomials in z/R with real coefficients so that 

in 

°‘ 15) H m (2 > * J 0 W E) V 2/E > 

which is the desired approximation* 

Remark 1 : The polynomials Q k are the imaginary analog of Chebyshev 

polynomials. They are orthogonal on the interval [— ± , i] with respect to the 
following inner product 

(3.16) <f,g> = -if f (w)g(w)( 1— | w | 2 ) ^ 

-i 

Remark 2 : It is apparent from (3.5) that in using the maximum norm we did 
not take into account the fact that the b k 's are decreasing. To do so 
requires us to consider the larger set of Gegenbauer polynomials. When the 
degree m is large it can be shown that the improvement thus achieved is 
negligible so that only for small values of m is the larger set relevant. A 
detailed analysis of the use of Gegenbauer polynomials will be carried out in 
a future paper which will deal with nonlinear problems. 
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4. Resolution in Time 

Let us define first the notion of resolution. The accuracy of a 
polynomial approximation is defined by its asymptotic rate of convergence as 
m (the degree of the polynomial) tends to infinity. Denote by [m o ,°°) the 
interval of the asymptotic behavior, m has then to be greater or equal to 
m Q in order to have resolution. This is a necessary condition but not 
sufficient. For example, if the relative error is of order 1, the results are 
meaningless and we have no resolution. Therefore, we define the condition of 
having a meaningful resolution as one in which m ^ m^ and the relative error 

norm is less than 10%. To be precise, assume that for m € [n^, 00 ), the 

minimal m which achieves this accuracy is m Q , we then say that a necessary 
and sufficient condition for resolution is 

(4.1) m>“ 0 

Applying the above definition to our case means that one has to apply the 

spatial operator tG^, m Q times in order to resolve N modes of the exact 
solution of (3.1) at time level t. 

Let us see what is m in the case of Chebyshev polynomials 

o 

approximation. 

Using the results from the previous section we have 

00 

(4.2) e z = k I Q Vk (R)Q k (z/R) 

It is known [10] that J k (R) converges to zero exponentially fast when k 
increases beyond R. It implies that the interval of asymptotic behavior is 
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[m o ,°°) while m 0 = [R] • Because of the exponentially rate of convergence 

m is close to m and when R is large we can consider m as nu for 
u w o o 

any practical use* Thus, we obtain that in order to resolve N modes one has 
to use the spatial operator at least | a t ( N— 1 ) | times. 

For comparison let us analyze the resolution qualities of the leap-frog 
scheme. This scheme is a typical explicit scheme which evaluates the 
numerical solution at the n + 1 time set, using data from the two previous 
time levels 

(4.3) V n+1 = V 11 " 1 + 2AtGV n 


A straightforward eigenvalues analysis implies that there are two solutions 
for the amplification factor of each mode w^ . 


(4.4) 


with 

(4.5) 


AtX, 

k 



= AtX, 
k 



2 


+ 1 


X, = iak 
k 


The scheme is stable when |akAt| 1 and we get 


(4.6) |„ 1>2 | = l 

which means that the error of the scheme is only a phase error. 

Let us assume that we choose the initial data at the first two levels in 
such a way that only y 1 is relevant. Therefore 
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(A. 7) 
or 

(A. 8) 
where 
(A. 9) 


t 1 = uV 


± t . 

V 11 = e in %° = e 


'I' = tg 


-1 


(i-e 2 ) 1/2 


= e + + o(e 5 ) 


is the phase shift of the numerical scheme after one time step. The quantity 
0 = akAt is the phase shift of the exact solution after t = At. Hence, the 
phase error at time level t is 


(A. 10) 


A0 = ^ ^ + O(0 5 ) = tk 6 a A t 2 + 0(At A ) 


The largest mode is w n (n = so that the maximum phase error is 


(A. 11) 


A0 = S(air)^ + 0(At A ) 
max 6 Ax 


This scheme is obviously second order in time and error E is 


(A. 12) 


E - |e“-e i(e « 6) | - |l-e“ e ||e 18 | . 


When A0 > ir, decreasing At would not necessarily decrease the error. Thus 

resolution is achieved when A0 is at least less than it . Therefore since 

max 

(A. 11) we obtain 


3 

a 



-j + 0(At 4 ) 


(A. 13) 
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or, using the notion of mQ defined at the beginning of the section, 


(4.14) 


_ it ( a» t -) 3/ 2 _ 1^ 

(6tt) 


“0 ~pj 2 ^Ax^ 


^(taN) 372 + 0(4t“) 


and the sufficient condition for 10% error results in 


(4.15) 


|A0| < 10 


-1 


Hence, using (4.10) we have 


m 0 = ^ 1/2 (taN) 3/2 + 0(At 4 ) 


Thus we have obtained the following result: in order to get a resolution of 

N modes by the leap-frog type scheme one has to operate with tG N at least 
^ \j 3/2 

(y) 2 | (taN) | times, a requirement much more stringent than the previous 

result of | ta(N-l) | operations for the Chevyshev approximation. For 
example, when t = 2tt , a - 1 and N = 32 one has to apply the spatial 
operator, in the leap-frog case, approximately 1300 times compared to 100 
times in the Chebyshev approximation case. 

As stated previously, for any practical use one can identify Chebyshev 
polynomials approximation as the "best polynomial". Thus we conclude that 
from resolution point of view the scheme based on these polynomials yields the 


best results. 



14 


5. Stability 


In the last section we discussed in detail the notion of resolution, 
is clear from the above discussion that resolution implies stability, 
fact, since 


(5.1) 



< const 


It 

In 


H^GnO must be bounded independently of m and N. 

The converse is not true in general. Consider for example the leap-frog time 
difference scheme. It has been shown in Section 4 that in order to get 
resolution we need 


m > (j) 1/2 (taN) 3/2 


or, equivalently, 

(5.2) 


aAt . ( 3 i 1/3 1 . 
At — ''5m'' ir ’ 


this is in contrast to the stability condition 


(5.3) 


a*At , N 1 _ _1 
St" — N-l IT TT 


for the leap-frog scheme [5] which allows a much larger At. 

When a time step At is chosen based on the condition (5.3) rather than 

(5.2), one may get meaningless results in spite of the stability of the 
scheme. To illustrate this we solved the equation 


u - u = 0 0 <x<2tt 

t x — — 

(5.4) 

Uq(x) = cos kx 1 k _< 7 

numerically. 
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The exact solution is 

( 5 * 5 ) u(x,t) = cos(k(x+t)). 

Using a grid of 16 points in space assures us that the error at time level 
t comes solely from the time discretization. 

The results are: 


K 

1 

2 

3 

4 

5 

6 

7 

L2ERR0R 

.9236 x 10 -2 

.767 x 10 -1 

.2722 

.6070 

1.291 

2.096 

.9453 


t = 3.625; Ax = .3927; At = Ax/tt = .1250 

These computations illustrate the above claim: meaningful results are 
achieved only for 1 < k < 3, while for 4 < k < 7 the results are 
meaningless despite the fact that we have used a spatial approximation that 
resolves all the modes exactly. 

The conclusion is obvious: for nonstiff problems, the important property 
is resolution rather than stability. It is inefficient to use a scheme which 
is stable but does not resolve all the modes. A scheme with less modes and 
the same degree polynomial in time will produce the same results. In the 
leap-frog case, for example, any results achieved by using the maximum time 
step allowed by the stability condition could be achieved with less amount of 
work by using coarser grid and the same time step. 


Condition (5.2) can also be written as 


a 



thus, it is obvious that in order to get resolution in the leap-frog case t 

has to be proportional to (Ax)^ and not to Ax as required by 

3/2 

stability* A same proportion between At and Ax is needed to get 

resolution for any scheme second order in time, and it does not matter if the 
scheme is stable or not for At proportional to Ax. 

For any scheme of order P the truncation error can be written as 


E = c 


At 


P+1 


N P+1 ♦ 


r\( a P+2...P+2\ 
0[At N J ; 


thus, the overall error in time t is 


t , P+1 P+1 _ r P+2..P+2^ 

E = ~ • c • At • N + O^At N J 


or 


E “ t • c • ir 


P+1 At 


Ax 


P+1 


+ 0 


At 
, Ax 


P+1 

P+2 


It follows that for a scheme of order P , At has to be proportional to 
P+1 

p 

Ax in order to get resolution. 

Considering the result of Lemma 2 and using the relations At = t/m, 

Ax = ir/n the requirement for resolution in the Chebyshev polynomial case is 
equivalent to 

a*At x N 1 1 
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which means that At has to be proportional to Ax. Thus, this new algorithm 
can be regarded as*a limit case of finite order schemes. 

We would like to mention here another important result which follows this 
stability discussion. The stability condition (5.3) for the leap-frog type 
schemes is much more stringent than the C.F.L. condition of a similar 
algorithm based on finite difference approximations despite the fact that 
spectral approximation uses all the previous mesh points. Thus, it was 
regarded as an artificial condition which may be overcome by properly designed 
time algorithm. Observing the fact that this stability restriction is exactly 
(5.6) which is the resolution condition for the orthogonal polynomials 
algorithm we conclude that this severe stability conditon is an essential one 
which can not be violated. The reason is due to the fact that the spectral 
radius of the operator is increasing. For example, in the leap-frog type 
algorithm the eigenvalues of Gjj, using finite difference approximation in 
space, are 


(5.7) 


. sin(KAx) 
k 1 Ax 


-N < k < N 


and the spectral radius is 


I . sin(KAx) | .IN 

r wn = max i x < — = — 

FD 1 Ax 1 Ax it 


On the other hand, for spectral (in space) approximation, the eigenvalues are 


(5.8) 


X, = ik 
k 


-N < k < N 


hence, the spectral radius is 
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(5.9) 
so that 


SP 


max | ik | = N - 1 
k 



N - 1 

TT 5 


ir 


6. Accuracy 

According to (4.4) we have the following expression for the coefficients 


\ * (1) Vk< E) c 0 ‘ 2 ' <Tc ' 1 k > 1 


Bessel functions satisfy the following inequality [10] 


(6.1) J (nut) I < 

m — 


<|)*exp 




l-t/l -* 2 


m 


kl < i 


Define 

( 6 . 2 ) 


4>*exp( 


l+/l-<J) 2 


(f»| < 1 


so that 



4> > 0 


* < 0 


where 


8 = /l-4> 2 
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Thus a is monotone decreasing for -1 < <j> < 0 and increases monotonically 
for 0 _< <j> < 1; and also 

(6.5) a(0) = 0, a(l) = a(-l) = 1 

Hence it follows that 

( 6 *6) 0 < a < 1. 

In our case m = R which implies 

(6.7) * = R / m = ta^i- ~ ta£ 

m m 

Thus, refinement of the approximation in space and time keeping the same 

proportion of N/m reproduces a scheme whose error in time converges to zero 
as 

(6.8) a m — — =► 0 0 < a < 1 

m ■+■ 00 — 

Hence, we have produced with a scheme which has spectral accuracy both in time 
and space. 

We would like to point out an interesting result which can be concluded 
from this analysis. When T is large m is large as well since it has to be 
greater than R in order to have resolution. According to (6.8), once we 
obey this requirement of resolution ( d> <1 ) the time error is negligible, and 
the error of the numerical solution of the (2.1) comes only from the spatial 


approximation. 
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7. The New Algorithm 

In this section we describe the actual construction of the algorithm. 
In order to obtain our scheme we use the expansion (4.2) 

m 

e Z - I c. J (R)Q, (z/R) 
k=0 

c 0 = 1, c k = 2 k > 1 • 

Substituting G N t for z results in 

(7.1) V(t) = e^V = l c k J k (R)(Q k (G N )v°) 

k=0 

where ~G = . Using the recurrence relation (3.14) and (3.11) we get 

N R N 


(7.2) 

MSP - v o • 

Using (7.2) in (7.1) enables us to compute V(t) by operating with G^ m 
times. It is obvious that because of the use of the recurrence relation this 
scheme may be regarded as a two level scheme. Therefore, it has the 
disadvantage of requiring extra memory. This disadvantage can't be overcome 
by converting (7.1) to a power series in G N and using Horner scheme to 
compute V(t) because huge roundoff errors result. 

A useful way to compute T(t) by a one level scheme is to calculate the 


roots of the polynomial 
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< 7 * 3 ) P(z) = I a J (R)Q (z). 

k=0 K K K 

Let us assume that the roots are 


(7.4) 


^1 *^2 * ' 




m 


Since c kJk(R) are real, every complex root appears with its conjugate. 
Rearranging (7.4) in such a way that the first 2p roots are p couples of a 
complex number and its conjugate we get 


(7.5) 


VV 


•*,u ,u ,U ,< 

P P 2p+l 


m 


Thus, (7.3) can be written as 


(7.6) 


P o rn 

p(z) = a n (l-a z+P.z J n (l-y z) 
i=l 1 i=2p+l 1 


while 


m/2 

= I 

k=0 


C 2k J 2k (R) 


g i 


1 



(7.7) 


a, 


2R e"i 




Hence we get 


(7.8) 


P 

p(z) = a n 
°i=l 


I 


«i(S) + 6 i(°n ) 2 


m 

n 

i=2p+l 


- 


v° 
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This is obviously a scheme that uses the minimal memory required for an 
explicit scheme. 

Our algorithm can be used as one step method by getting the solution at 
the final time t directly from the initial data. It can also be used as a 
marching scheme if one is interested in intermediate results. The size of the 
time step At depends only on the information one wants to get out of the 
numerical procedure. At enters instead of t in the expressions above and 
the parameters R,m are determined accordingly. In any case, the refinement 
of the algorithm is done by increasing the degree of the polynomial and not by 
decreasing the size of the time step. 


8* Variable Coefficients 

In the variable coefficient case, the operator G is 

(8.1) G = a(x)|~^ 

It is approximated in the numerical procedure by the matrix G^ which is a 
multiplication of two matrices 

(8.2) S ' "n • °N 

where is a diagonal matrix whose elements are 

^Vij = a ( X j) 6 ij 


(8.3) 
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and Djj is a skew hermition matrix which approximates the derivative 
spectrally. It is clear from this representation that G N is no longer 
normal. Nevertheless the main results of our approach are still valid. In 
the non-normal case, the set of eigenvectors is not always complete and one 
cannot use (3.4) - (3.5). 

in order to justify our approach in the non- normal case we use the 
following definition for a function fo a matrix A [3]: 


(8.4) 


f (A) = l 
k=l 




z. .+f (1) (X. ) 


k'' Z k2 


+ *** +f ^kKn^ 


where X fc are the eigenvalues of the matrix A and m^ is the multiplicity 
of X^ in the minimal polynomial of A. The matrices z^j are completely 
determined when A is given and do not depend on the choice of the function 
f. Expression (8.4) has a meaning when f and its required derivatives are 
defined on the spectrum of A. In our case f is the exponent function which 
is a well defined analytic function. Hence, using definition (8.4) it is also 
obvious that in the general case approximating the exponent matrix is 
equivalent to approximating the scalar exponent and its derivatives on a 
domain which includes all the eigenvalues of the matrix. 

In the constant coefficients case the domain is 


(8.5) I = [-iaN.iaN] . 

The following theorem implies that this result is valid also in the variable 
coefficients case when a(x) doesn't change sign in the interval and 
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(8.6) a = max|a(x)| x€[0,2tt] 

x 

Theorem : If a(x) > 0 and X^ is an eigenvalue of A N D N then X^ € 

(I is defined by (8.5), (8.6).) 

Proof : Define the following inner product 


(8.7) 



Then we have 

( 8 . 8 ) 


[ 1 v( u n) J.-i [“n’WJ -l 


where (u^,D u ) is real; thus 


(8.9) 
Hence 

( 8 . 10 ) 


& 

( U N ,D N U n) = ^ d n u n ,u n^ = ^ U N’ D N U J = _ ^ u n ,d n u n^ 

(vvJ = 0 


Using this result in (8.8) we get 


( 8 . 11 ) 



thus 

( 8 . 12 ) 



= const, 


Assuming that w^ is the eigenvector of A^D^ corresponding to X^ we 


can 


use it as the initial data so that 
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(8.13) 


U N 


= e 


w 


w k = e 




w,. 


hence 

(8.14) 


l- U N ,U N-l -1 


- [ 


V 


w k ,e 


V 


\] -1 
V 


2 (VkK 

e K 




and according to (8.12) this is constant. This implies that 


(8.15) R X = 0 

e k 

hence X^ is pure imaginery. In addition 


(8.16) 


max ' x kl < 'W < lA^I ID I = a*N 
lc 


and the proof is concluded. 

The proof is essentially the same when a(x) K 0 in the interval. 

The case when a(x) changes sign is more complicated and not much is yet 
known. [6] gives a proof of stability for the simple case when a(x) is a 
trigonometric polynomial of order 1. It implies that in this simple situation 


(8.17) -a < R X. < a a>0 

— e k — 

while a doesn t depend on N. Because of (8.16) there exists an ellipse 
whose larger axis is [-iaN,iaN] and the small one is [- 0 , 0 ] which contains 
the eigenvalues of A^D^. The theory of Chebyshev polynomial approximations 
guaranties convergence in this domain [9]. 
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Remark Numerical experiments lead us to the assumption that (8.17) is 
also valid in the general case that a(x) is any periodic function. 


9. Numerical Results 

In order to illustrate the spectral convergence of our scheme we shall 
consider the following scalar problem 

(9.1) u - a(x)u =0 0 ^ x £ 2n . 

C X 

In Table 1 we take 

(9.2) a(x) = g— — , u (x) = sin(2x+sin x) 

l + COS X o 

The exact solution to this problem is 

(9.3) u(x,t) = sin(2x+sin x+t) 

The numerical solution is computed at time leve T = 6.283 
N = Number of mesh points in space. 

M = The degree of the ploynomial approximation. 

The ratio of the L 2 errors illustrates very clearly the spectral convergence 


of the scheme 
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Table 1 


N 

M 

L»2 Errors 

Ratio 

8 

36 

1.605 x 10 -1 

9.2 x 10 3 

16 

72 

1.740 x 10“ 5 

4.6 x 10 7 

32 

144 

3.756 x 10" 13 


In Table 2 

we take 



(9.4) 

a(x 

:) « sin(x), u (x) = sin(x) 

o 


(a(x) changes 

sign in the 

interval) • 


The exact solution is 



(9.5) 


u(x,t) = sin( 2tg ^(e^gj)) 


The sollution 

is computed 

at T = 1.571. 




Table 2 


N 

M 

L2 Errors 

Ratio 

16 

18 

5.968 x 10 -2 

2.9 x 10 1 

32 

36 

2.031 x 10 -3 

8.7 x 10 2 


2.345 x 10 


64 


72 
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Table 3 inllustrates the resolution properties of the scheme. In this 
case we 

a (x) = - , T = 50.27, N = 16 

2 + cos x 

Because of the high resolution in space only time errors occur. According to 
the theory developed previously (Lemma 2) the degree of the polynomial 
approximation has to be at least R. We have 

R = at(N-l) = 50.27 x 15 = 754.05 


Table 3 


M 

L2 Error 


740 

1.120 



750 

5.981 

X 

10 -1 

760 

1.354 

X 

10 -1 

770 

1.476 

X 

io -2 

780 

1.048 

X 

CO 

1 

o 

840 

5.391 

X 

io- 13 


The result for M = 840 shows that, while the minimal M required for 
resolution is 755, increasing M by only 11% gives machine accuracy. 

In Tables 4 and 5 we compare our scheme to the leap-frog scheme. The 
model problem is (9.1), (9.2) with N = 32. In Table 4 we compute the 

numerical solution at different time levels. The results clearly shows the 
high accuracy of the Chebyshev polynomials approach. Another interesting 
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phenomenon which is illustrated by the table is the following: while the 

leap-frog case the error increased due to accumulation this is not so in our 
case. In contrast it decreased geometrically. The explanation is obvious 
using (6.7), (6.8). Increasing t and m, keeping the same proportion 

between them, results in <f> = constant and therefore (6.8) is valid. 

In the following table 
C.P. = Chebyshev polynomials. 

L.F. = leap-frog. 


Table 4 


Time 

M 

L2 Error 

(L.F.) 

L2 Error 

(C.P.) 

1.571 

35 

5.247 x 

10 -4 

8.726 x 

10" 5 

3.142 

70 

1.108 x 

10 -3 

2.084 x 

10 -7 

6.283 

140 

2.184 x 

IQ" 2 

1.429 x 

10" 13 

In Table 5 we 

compare our scheme to the 

leap-frog scheme from 

the point of 

view of the amount of work needed to achieve a certain degree 

of accuracy. 

The L2 Error is computed at time level T = 

6.283. 





Table 5 




L2 Error 

u(L. F. ) 


u(C.P. ) 


1.0 : 

x 10" 4 

580 


no 


1.0 : 

x 10~ 6 

3480 


117 



-8 





1.0 : 

X 10 

56000 


122 
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